
# read data and load the GAPIT function
# 这个GAPIT软件有三个主要功能，第一个是群体遗传学分析（我们叫他 population anlysis），第二个是全基因组关联分析（GWAS），第三个是全基因组选择（GS）
# 这三个部分分成三个软件，各自的运行窗口相互独立，但是设计风格上保持相似。
# 总体上的功能是每个软件都需要窗口式操作，可以用鼠标设置参数，选择模型，设置输出格式等。
# 
rm(list=ls())

#每个软件都需要进行的步骤，以下是引用GAPIT数据包和所依赖的包。
library('MASS') # required for ginv
library(multtest)
library(gplots)
library(compiler) #required for cmpfun
library("scatterplot3d")

source("D:\\RProg\\emma.txt")
source("D:\\RProg\\gapit_functions.txt")


# 下面是加载数据，三个软件的加入风格可以保持相似。
# 首先是表型数据，这个可以不是必须的， , head = TRUE)
myY=read.table("D:\\RProg\\mdp_traits.txt",head=T)
# 下面的是hapmap基因型数据，有他这个一个了，就不用后面的两个了
myG=read.table(file="D:\\RProg\\mdp_genotype_test.hmp.txt", head = FALSE)
# 下面是numeric格式的基因型数据，分为两个文件，这两个文件相互依存，有一必须有二，有二必须有一。
myGD=read.table("D:\\RProg\\mdp_numeric.txt",head=T)
myGM=read.table("D:\\RProg\\mdp_SNP_information.txt",head=T)


########################### PA ###########################
# PA 的第一个功能是实现数据转换和群体分析
# 如果直接放入G（hapmap）的话，那就可以转换成myGD和myGM（numeric）格式
myGAPIT=GAPIT(
#G=myG,   #如果有这个了，就不用后面两个数据
GD=myGD,
GM=myGM,
kinship.algorithm="VanRaden",# 这是用来计算 分子亲缘关系矩阵的参数，里面可以是VanRaden，Loiselle，Zhang,Separation，默认是VanRaden。
PCA.total=3, #这是告诉模型要计算的PCA数量，这里可以是0以上，小于个体数量的任何证书
PCA.3d=TRUE, #这是做3D PCA的参数，值为TRUE时，开始做，默认是FALSE
NJtree.group=4 #这个是做NJtree分析的关键数值，默认为NULL，数值设置多少，就是将群体聚类成几个组
) 


########################### GWAS ###########################
setwd('D:\\RProg\\22')
myGAPIT<- GAPIT(
Y=myY[,c(1,2)],
#GD=myGD,#Genotype
#GM=myGM,#Map information
G=myG, #GWAS 和PA 都支持两个基因型文件（GD,GM）输入，或者只放入一个G
PCA.total=3, #设置PCA分析中使用的个数，最小（默认）是0
Inter.Plot=TRUE, #选择是否进行交互分析
Multiple_analysis=TRUE, #选择是否进行多重分析比较
group.by=50, #设置压缩比例，
#Random.model=T, #选择是否使用随机模型进行后续验证
model=c("GLM","MLM","CMLM","SUPER","MLMM","FarmCPU","Blink"),#GWAS中共有7种模型，这里可以做成选择对应好代码就行
file.output=T #选择是否输出结果
)



########################### GS ###########################

set.seed(198521)
nfold=5 # 选择交叉验证的比例
sets=sample(cut(1:nrow(myY ),nfold,labels=FALSE),nrow(myY ))

j=1
training=myY[,c(1,2)]
training[sets==j,2]=NA
training_index=is.na(training[,2])
testing=myY[training_index,c(1,2)]

myKI=GAPIT.kinship.VanRaden(snps=as.matrix(myGD[,-1]))
# Time=timer(time=Time,infor="kinship.end")
kinship=cbind(as.data.frame(myGD[,1]),myKI)
colnames(kinship)=c("Taxa",as.character(myGD[,1]))

myGAPIT<- GAPIT(
Y=myY[,c(1,2)],
#GD=myGD,#Genotype
#GM=myGM,#Map information
G=myG,
PCA.total=3,
#KI=kinship, # GS中，可以放入基因型文件（2个,GD,GM）或者hapmap格式文件（1个,G），或者KI，都可以直接进行GS。
Inter.Plot=TRUE, #这里也涉及到GS中的交互分析结果
group.by=50,
#Random.model=T,
# model=c("gBLUP","cBLUP","sBLUP"), #这里也是选择GS中方法，可以做成选择
model=c("gBLUP"),

# testY=testing, #这里可以做个判断，选择是否放入
file.output=T  #一样判断是否输出结果
)


